{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ad2687b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3f097ec",
   "metadata": {},
   "source": [
    "#初始化数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a3b1866c",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng()\n",
    "\n",
    "xn = np.random.randn(100,1)\n",
    "\n",
    "x = signal.lfilter([1,1], 1, xn)\n",
    "x = np.array(x).flatten()\n",
    "y = signal.lfilter([1,-1], 1, xn)\n",
    "y = np.array(y).flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "765c9cce",
   "metadata": {},
   "source": [
    "#定义函数detrend_none\n",
    "\"\"\"\n",
    "    Return x: no detrending.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : any object\n",
    "        An object containing the data\n",
    "\n",
    "    axis : int\n",
    "        This parameter is ignored.\n",
    "        It is included for compatibility with detrend_mean\n",
    "\n",
    "    See Also\n",
    "    --------\n",
    "    detrend_mean : Another detrend algorithm.\n",
    "    detrend_linear : Another detrend algorithm.\n",
    "    detrend : A wrapper around all the detrend algorithms.\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "80d58f88",
   "metadata": {},
   "outputs": [],
   "source": [
    "def detrend_none(x, axis=None):\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d182899",
   "metadata": {},
   "source": [
    "#定义函数M_xcorr\n",
    "r\"\"\"\n",
    "    Plot the cross correlation between *x* and *y*.\n",
    "\n",
    "    The correlation with lag k is defined as\n",
    "    :math:`\\sum_n x[n+k] \\cdot y^*[n]`, where :math:`y^*` is the complex\n",
    "    conjugate of :math:`y`.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x, y : array-like of length n\n",
    "\n",
    "    maxlags : int, default: 10\n",
    "        Number of lags to show. If None, will return all ``2 * len(x) - 1``\n",
    "        lags.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    lags : array (length ``2*maxlags+1``)\n",
    "        The lag vector.\n",
    "    c : array  (length ``2*maxlags+1``)\n",
    "        The auto correlation vector.\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "214fc069",
   "metadata": {},
   "outputs": [],
   "source": [
    "def M_xcorr( x, y, normed=True, detrend=detrend_none,\n",
    "           maxlags=10, **kwargs):\n",
    "    \n",
    "    Nx = len(x)\n",
    "    if Nx != len(y):\n",
    "        raise ValueError('x and y must be equal length')\n",
    "\n",
    "    x = detrend(np.asarray(x))\n",
    "    y = detrend(np.asarray(y))\n",
    "\n",
    "    correls = np.correlate(x, y, mode=\"full\")\n",
    "\n",
    "    if normed:\n",
    "        correls /= np.sqrt(np.dot(x, x) * np.dot(y, y))\n",
    "\n",
    "    if maxlags is None:\n",
    "        maxlags = Nx - 1\n",
    "\n",
    "    if maxlags >= Nx or maxlags < 1:\n",
    "        raise ValueError('maxlags must be None or strictly '\n",
    "                         'positive < %d' % Nx)\n",
    "\n",
    "    lags = np.arange(-maxlags, maxlags + 1)\n",
    "    correls = correls[Nx - 1 - maxlags:Nx + maxlags]\n",
    "\n",
    "\n",
    "    return  correls, lags\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b944acc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArrElEQVR4nO3de3hV9Zno8e9LJIJCAS+NGFCxcgsgCUlBgSKhKOjpCGod8bRWWjzp2OPTU9tBoTrCY2tFmbYznV6mjsOonQ60AkYsWKpA2lK1gnJJICARHCUg18QkmkAI7/ljr+zZe2dfs/YtWe/nefaTtX7rst/122uvN+u31l4/UVWMMcZ4V49MB2CMMSazLBEYY4zHWSIwxhiPs0RgjDEeZ4nAGGM87pxMB9AZF110kV5xxRWdWvbjjz/m/PPPT25ASWBxJcbiSozFlZjuGtdbb711XFUv7jBBVbvcq7i4WDtr06ZNnV42lSyuxFhcibG4EtNd4wK2aphjqjUNGWOMx1kiMMYYj7NEYIwxHtclLxaH09raysGDB2lpaYk6X79+/aiurk5TVPGzuBLTFePq1asXgwYNomfPnmmOypjouk0iOHjwIH379uWKK65ARCLO19jYSN++fdMYWXwsrsR0tbhUlRMnTnDw4EGGDBmSgciMiSwpTUMiskxEjopIVYTpIiI/EZEaEdkpIuMCpt0tIvuc192djaGlpYULL7wwahIwJlXqPjnNnsMNHPiojT2HG6j75HTQ9PrmVo6d7sm+wyeZtGQj5dtqMxSpMR0l6xrBM8DMKNNvBIY6rzLgFwAicgGwCJgAjAcWiciAzgZhScBkQt0np6mta+Z021kATredpbau2Z8M2qe3nlUEoba+mYWrKy0ZmKyRlESgqn8CTkaZZRbwnHMr6xtAfxEZCMwAXlHVk6paB7xC9IRiTNY58lELZ0Me535WlSMftUSc3tzaxtL1e9MWozHRpOsaQT7wQcD4QacsUnkHIlKG72yCvLw8Kioqgqb369ePxsbGmIG0tbXFNV+6WVyJyaa42s8EwpU3NjZGnF5b39xhP06VpqamtL1XIiyuxKQqri5zsVhVnwKeAigpKdGpU6cGTa+uro7r4mFXu8iYaRZXbLlNDWEP9rk5Pejbt2/E6fn9exO6H6dKRUVF2t4rERZXYlIVV7p+R1ALDA4YH+SURSpPufJttUxaspEhC9Ym9eLd4cOHmTNnDiUlJQwbNozS0tKkrDeSDRs28OUvfzml79HVxVNHBw8e5De/+Y1//LXXXuORRx6Ja/15/XrRI+T6VA8R8vr1iji9d88c5s8YHtf6jUm1dCWCNcBXnLuHrgE+UtXDwHrgBhEZ4FwkvsEpS6nybbUsXF1JbX0zCkm9eHfXXXdxyy23sHXrVt555x1+8pOfuA84ih07dlBYWOhqHW1tbckJJo1CY462DfHU0YYNG3j77bf94xMnTuTRRx+NK5YB5+WSP6C3/2aF3Jwe5A/ozYDzcoOm57adQVDy+/fm8VvHMLsobCuoMWmXrNtHlwOvA8NF5KCIzBORvxORv3NmWQfsB2qAfwO+AaCqJ4HvAVuc16NOWUotXb+X5tbgA0cyLt61tbVRUVHBdddd5y8bM2YMACtXruSaa65h7NixTJ48mWPHjgFw++23c99993HDDTdw+eWXs3nzZu666y6GDRvGvHnz/Ou58847ueOOOxg/fjyXX345a9euBXwHubFjxwJw4MABZs2aRUlJCePHj2fvXt/2HDp0iNtuu42ioiJGjBjBm2++ye23387Xv/51rrnmGh5//HEA9uzZw7Rp0ygsLGT69OkcP34cgGeffZbi4mKuvvpqJk+eHLEsUeHiihZHYMw//OEPO2xDpO0PrKNwn8PmzZv59re/zcqVKyksLGT//v3cfvvt/PnPf44az6233srDDz/MlClTGDviKra9/id65QgjBn7KnwTaDTgvlxHH3iO/4Rh/WTDNkoDJLuGeRJftr3BPH929e3dcT99raGjQKx78nV4e5nXFg7+Lax3RzJgxQz/96U9rWVmZbt682V9+/Phx//DixYv1pz/9qaqqDh8+XH/4wx9qQ0ODPvbYYzps2DA9dOiQtra2al5enra0tKiq6siRI3XBggWqqvrnP/9ZP/vZz6qq6tixY/Xo0aN6+vRpnTZtmtbU1Kiq6tq1a3Xu3Lna2tqqV199tb700kuqqvrxxx9rQ0ODDh8+XP/hH/7BH1NLS4sWFBTotm3bVFV1yZIl+t3vfldra2t15MiReurUKVVVraur04aGhg5liYoUV6Q42uuqPebQbYi0/YF1FO1zmDFjhlZWVvqnjRgxQuvr66PGc9VVV+nSpUtVVXX16tV625wv6d5D9ZE3essW3f3yywnXVTJ016dppkp3jQt7+uj/uLR/74TKE/Hyyy+zatUq+vXrx8yZMykvLwfgmWeeYfz48YwdO5af//zn9OrVi5aWFurr6/nWt74F+H4HMW/ePAYOHMg555xDTk4Oubm5tLS0cOzYMRYtWgRAQUEBdXV1tLa28tFHH3HxxRdTXl7Orl27uO222ygsLOSBBx6gV69elJeXM3LkSL7whS8AcN5559GzZ09OnjwZ1AZeXl7O5MmT/U0oBQUFHD16lJycHJqbm/nOd77D1q1b6d+/f9iyQNOnT2f06NEdXi+++GLQ+4XG1bdv34hxtLS0BMUcOh5p+wPrKNLnALB3715GjBjhX/fp06fp169fxHg++eQTPvroI+6//37A94iTvv2C68GYrqLL3DWUTPNnDGfh6sqg5qFkXbwTESZPnszkyZOpq6tj586dNDQ08Oabb7Jx40b69OnDlClTGDVqFLt27WLcuHH06OHLxzt27ODee+8FfBcvL730UkSEqqoqhg4d6j9ovf3224wdO5bq6mpGjhzpX/axxx4Lak4CePjhh7nmmmuCynbt2sWECRM455z/+fh3797tb8YCqKyspKCggPPOO4+qqipeeuklysrKuOeee/jGN74Rtqzdq6++GrOetm/f3iGuaHGExlxdXR00Hmn7d+7c6a+j5557LuzncPz4cfr16+df165duygoKIgaz+7duykuLiYnJ8f/PsNGjIy53cZko26ZCD78wQ84Vb0n7LQzbW0U5eTwnd6X8fPzR3Mytw95bZ/wf07upOhH/8V/R1jnuSNHcMl3vxv1fdevX09paSm5ubkcPXqUzZs3s2zZMlavXs3EiRPp06cPq1at4rXXXmPMmDE8//zz/rZr8B1Mrr76asB3YAscfv/992lpaaGtrY1Fixbx5JNPBrV9Dxw4kPXr1/PVr36VHj16UFlZyejRo7nkkkvYsWOH/z2OHTtGZWWlf93t8vPz2b59OwD79+/nV7/6FZs3b6ampoaioiLmzJnD7t27aWlpYd++fQwdOjSoLFHh4rr44osjxvHSSy8Fxbx79+6g8UjbH1hHlZWVYT+H6upqLr30Uv+6AusnUjxr1qwJugC9c+dOvnrf9QnXgzHZwJNNQwDXN7/PC8fX8cdDv+W3R37H9c3vu17nypUrGTlyJGPHjuULX/gC3/ve97j22muZO3cuP//5zxk/fjzbtm3jyiuv5Pzzz6eystJ/MGlpaaG5uZkBA3xP2AhNCrfeeisTJkzgs5/9LPfeey+TJk0Kuhvma1/7GmfPnmXkyJEUFhbyxBNPICLMnTuXI0eOMGrUKAoLC3n99dfDJoK77rqLQ4cOMWbMGObMmcOyZcu48MIL+cd//EeGDx/OuHHjOHDgAN/4xjd47LHHOpQlKlxc0eIIjXnXrl1B45G2P7COIn0OI0aM4Pjx44wePZrXXnst6L2ixROYCKqqqhg2oiDhejAmK4S7cJDtL7cXi7NRtLimTJmie/bsSWM0/6Mr1lem1BxttIvFCbK4EmMXiz3s3XffZejQoZkOwxjTTXXLawTdzcGDBzMdgjGmG7MzAmOM8ThLBMYY43GWCIwxxuMsERhjjMdZIjDGGI+zRGCMMR5nicAYYzzOEoExxnicZxNBKrqqfPrppyksLKSwsJAePXr4h++//342bNjAXXfdlYTI4+Om68VQbrrDDI0j01LdbaUxXVGyeiibKSJ7RaRGRBaEmf5jEdnuvN4RkfqAaW0B09YkI55YUtVV5T333MP27dtZu3YtgwcPZvv27Wzfvp0f//jH7Nixg6KiouRsgCNa94xuul4M5aY7zNA4UiGbuq00pitynQhEJAf4GXAjUADcKSJBj2FU1ftVtVBVC4F/AVYHTG5un6aqN7uNJx6p6qqyXVVVVdAz7MF3APrwww+ZMmUKl112mf+Z/e3dK1533XVB3SvG011jpO4ZY3W9GKmLyEjdacbq6hHg2muv5cCBAwDU1tZSXFwcNo5I25Vot5Vf+cpXguohG7qt/NHjj/KlW24M+nyN6RLCPYkukRdwLbA+YHwhsDDK/K8B1weMNyX6ntncVaWq6tKlS/XBBx8MKhs7dqw++eSTqurr1nDu3LlB3Ss2NDT4u1eMt7vGaN0zRup6MVIXkarhu3FsaGiI2dVjW1ubDhw4UM+ePauqquvWrQsbR6Tt6ky3lUOHDg3qajMbuq1csOj7uvdQvf/z7cCePtqBxZWYVD19NBkPncsHPggYPwhMCDejiFwODAE2BhT3EpGtwBlgiaqWR1i2DCgDyMvLo6KiImh6v379aGxsjBlsW1sbl3zqXA43nOow7ZJPnRvXOmLZtm0bpaWl/nW1trZy/PhxysrKaGxspKGhgfPOO4/ly5dTVVXFLbfcgqrS1tbGxIkTWb58ORMmTOAzn/kMjY2NDBkyhHXr1nHs2DFOnDjB/fffT2NjIy+88IJ/eYAzZ84wceJEGhsb2bNnD/n5+TQ2NtLS0kJLSws9evRg+fLlXHXVVVx33XVB29rY2Mgvf/lLVq1axenTpzly5AiLFi3yd6fZq1eviPPs2LGDyy67jKamJgC2bNnCsGHDOsSxatWqsNsVKaZo9VBXV+evh5aWlrjq5eTJkzG3JVbdhYvnyJEj1NfX85V77kVR/+cbui/1df6G7rvp0NTUlJH3jcXiSkyq4kr300fnACtVNbBd5nJVrRWRK4GNIlKpqu+GLqiqTwFPAZSUlOjUqVODpldXV9O3b9/QxTpobGzkwRtHhu2q8sEbR8a1jlj27NnDAw884F/Xzp07KSwspF+/fgDs27ePcePG8c477/CDH/yAefPm0djY6J9/0aJFjBs3zj9eU1PD2LFjef/997nmmmv8ndcELh/o+PHj9O/fP2i+0aNH07dvX/bu3cvnPve5Dtv53HPPsWPHDv74xz/6u3EsLi6mpqaGgoIC+vbtG3GeAwcOUFhY6F9nZWUlZWVlnDp1KiiOAwcOhN2uSDFFmv/999+npKQkaPviqZedO3fG3JbQmAPrLlI8H3zwASUlJeSem0vbmTb/5xtpXwrdd9OhoqIiI+8bi8WVmFTFlYxEUAsMDhgf5JSFMwf4v4EFqlrr/N0vIhVAEdAhESQi3q4q/63XCI727s+nk9RVJcDZs2fZt2+fv59cCG6XBt8BadasWbS0tPi7VwT83SvG211jpO4Z33vvvYhdL0bqIjJSN46//vWvY3b1WFVV5e/Avrq6mrVr1/LTn/60QxyRtmv58uUJd1s5evTosNsXrV7S2W3lrFmzOuwbxmSrZNw1tAUYKiJDRCQX38G+w90/IjICGAC8HlA2QETOdYYvAiYBu5MQU0zXN7/Pb+v+QEUSu6oE33+JgwYNIjc3118W2P8w+C4mjx49Oqh7xUmTJvm7V4y3u8ZI3TNG63oxUheRkbpxrKqqitnV44wZM/j973/Pl770JZ5//nkuvPBC8vLyOsQRabs6023lqFGj/PUQb72ks9vKwERlTNYLd+Eg0RdwE/AOvv/kH3LKHgVuDphnMb5rAIHLTQQqgR3O33nxvJ/XuqrMJIsrftZVZeIsrsRk88ViVHUdsC6k7JGQ8cVhlnsNGBNabowxJn08+8tiY4wxPt0qEfjOfIzJTgpw9mymwzCmg26TCHr16sWJEycsGZispKqcOHOGXjU1mQ7FmA7S/TuClBk0aBAHDx70P/IgkpaWFnr16pWmqOJncSUmG+M61niKs2fP0lbfO+z0Xm++yaDFi+Gb30xvYMbE0G0SQc+ePRkyZEjM+SoqKpL+8LdksLgSk41xLf7l69TX17P+wRvDz1BQEL7cmAzrNk1DxhhjOscSgTHGeJwlAmOM8ThLBMYY43GWCIwxxuMsERhjjMdZIjDGGI+zRGCMMR5nicAYYzzOEoExxnicJQJjjPG4pCQCEZkpIntFpEZEFoSZPldEjonIdud1T8C0u0Vkn/O6OxnxGGOMiZ/rh86JSA7wM+B64CCwRUTWqGpo38O/UdX7Qpa9AFgElOB7XPtbzrJ1buMyxhgTn2ScEYwHalR1v6qeBlYAs+JcdgbwiqqedA7+rwAzkxCTMcaYOCXjMdT5wAcB4weBCWHmu01EpuDr5P5+Vf0gwrL54d5ERMqAMoC8vDwqKio6FWxTU1Onl00liysx2RhXfX0zbW1tEeOa6vzNRNzZWF9gcSUqVXGlqz+Cl4DlqnpKRL4OPAtMS2QFqvoU8BRASUmJTp06tVOBVFRU0NllU8niSkw2xvWLvb7+CGLFlYm4s7G+wOJKVKriSkbTUC0wOGB8kFPmp6onVPWUM/o0UBzvssYYY1IrGYlgCzBURIaISC4wB1gTOIOIDAwYvRmodobXAzeIyAARGQDc4JQZY4xJE9dNQ6p6RkTuw3cAzwGWqeouEXkU2Kqqa4BvisjNwBngJDDXWfakiHwPXzIBeFRVT7qNyRhjTPySco1AVdcB60LKHgkYXggsjLDsMmBZMuIwxhiTOPtlsTHGeJwlAmOM8ThLBMYY43GWCIwxxuMsERhjjMdZIjDGGI+zRGCMMR5nicAYYzzOEoExxnicJQJjjPE4SwTGGONxlgiMMcbjLBEYY4zHWSIwxhiPs0RgjDEel5REICIzRWSviNSIyIIw078tIrtFZKeIbBCRywOmtYnIdue1JnRZY4wxqeW6YxoRyQF+BlwPHAS2iMgaVd0dMNs2oERVPxGRe4EngTucac2qWug2DmOMMZ2TjDOC8UCNqu5X1dPACmBW4AyquklVP3FG38DXSb0xxpgsIKrqbgUiXwRmquo9zvhdwARVvS/C/D8FPlTV7zvjZ4Dt+PozXqKq5RGWKwPKAPLy8opXrFjRqXibmpro06dPp5ZNpVTE9dqhVla908qJFuXCXsJtw3oy8dKeGY8rGbIxrsf/2kxbWxsPTwwf19TSUgAqNm1KZ1hAdtYXdD6uZOzbqYgr1dzGVVpa+paqloSWJ6XP4niJyJeBEuC6gOLLVbVWRK4ENopIpaq+G7qsqj4FPAVQUlKiU6dO7VQMFRUVdHbZVEp2XOXbavnVhkqaW32J/kSL8qvqNgpGFjC7KD9jcSVLNsb1i72vU19fHzOuTMSdjfUFnYsrWft2suNKh1TFlYymoVpgcMD4IKcsiIhMBx4CblbVU+3lqlrr/N0PVABFSYipg/JttUxaspG5v/+YSUs2Ur6tQ4jdytL1e2lubQsqa25tY+n6vRmKyJjk8OK+nerjVzISwRZgqIgMEZFcYA4QdPePiBQBv8SXBI4GlA8QkXOd4YuASUDgReakKN9Wy8LVldTWNwNQW9/MwtWV3ToZHHK2Nd5yY7oKr+3b6Th+uU4EqnoGuA9YD1QDv1XVXSLyqIjc7My2FOgDPB9ym+hIYKuI7AA24btGkPRE4MX/IC7t3zuhcmO6Cq/t2+k4fiXlGoGqrgPWhZQ9EjA8PcJyrwFjkhFDNF77DwJg/ozhLFxdGbQD9e6Zw/wZwzMYlTHueW3fTsfxyxO/LPbafxAAs4vyefzWMeTm+D7i/P69efzWMUm7mGZMpnht307H8csTiWD+jOH07pkTVNad/4NoN7son6LL+jNhyAX8ZcG0bvtFMd7jpX07HcevtN4+mintO8kDK3dyuu0s+f17M3/G8G698xhjuod0HL88kQjAV5nL33yf+vp61j84LdPhGGNM3FJ9/PJE05AxxpjILBEYY4zHeaZpyHRUvq2Wpev3cqi+mUvtuonJIrZvppclAo9q/7Vi+73Y7b9WBOwLZzLK9s30s6Yhj/Lir61N12D7ZvpZIvAoL/7a2q32B38NWbA24Qd/Bc7rhYceumH7ZvpZ05BHXdq/t/8hVqHl8UhlG242tg+7aa5oX3a2M55NTR3Z+Dm63TdN4uyMwKPc/Fox8GmISuJPQ4z2n7XbdaeKm+aKbG3qiFXXbs+AOvs5evVJAJlkicCj3Dyvxc2BLdYBIlsPmm6aK7K1qSNaXbtNyG4+R689SygbWNNQHNycPmdjM0e79l8rAvzm69fGvZybA1u0A8TsovysPWi6aa7I1qaOaHUd63Nys+54dHbfTLXueiywM4IY3PxnlK3NHG65eRpirANEtj4p1k1zRbY2dUSra7cH8mz9HN3ozscCSwSOSO2h3bFt2C03B7ZYB4hsPWi6aa5oX7ZdtjR1RKtrtwfybP0c3Uj1scDNNRm3ktI0JCIzgX8GcoCnVXVJyPRzgeeAYuAEcIeqvudMWwjMA9qAb6rq+mTEFKp8Wy3b3q/ndNtZJi3ZGHRaVr6tlvnP76D1rK8z7Nr6ZuY/vwOI7xQ30ilfIsvW1jeT/8bGDqeL0U4nY51qRtvmWNMfLq9k+V8/oE2VnPXruHPCYL4/23cwm12Uz9b/Psl/vuE7dc8R4bbi/PDLigQtO3/G8KC6BujZQ/wHiFjrTmV9RYu7PbZIzRWx6jpw+C8Lgh8aFut94/mcY21zuPqKVdfRPqdYcSe0j4TsX/HUZ6zvcyrqK5XHgmjHoNlF+THrwy3XiUBEcoCfAdcDB4EtIrImpMvJeUCdql4lInOAJ4A7RKQAXx/Ho4BLgVdFZJiqBqdOl9pPy063nQU63r63eM2uoB0eoPWssnjNLvr17kl9c2uHdfbr3TNo3eFuK3SzbPuHH2k6ENeykbY52vTALzBAm6p//Puzx1C+rZZVb9UGTV/1Vi0ll18Qc1kAJKRCAsajrTuebZ6/cgetbQFfppU74qqvh8srY8cdQay6jibW+8azj3R2m2PVdbTPKZ64O7uPuNl349nmzn6nUnksiHYMao+rM/tXvERVY88VbQUi1wKLVXWGM74QQFUfD5hnvTPP6yJyDvAhcDGwIHDewPmiv2eJwlZXcRtjjPfIW6paElqajGsE+cAHAeMHnbKw8zid3X8EXBjnsgCISJmIbBURywDGGJNEXeb2UVV9CngKoKSkRLcmkA4mLdkY9va9/P69+cuCaRQ9+gfqPul42jbgvJ6cl3tO1GWjrRvo9LKxpn/4UQttYc7mckR49/GbMrbuVMYF0evzigVrO0wLnCfasp9ZuC5q3H7itI8EzBtrm6KJ9b6x1u1mm93UtZu4s3XfjbXNqdx33RyDEiGhzX2OZJwR1AKDA8YHOWVh53Gahvrhu2gcz7KuxbqDYdHfjKJnTnAN9cwRFv3NqJjLRpsea9nSEReHjbe9PNry4XZ2wF/uJu47JwwmnPbyVC0ba3qsZfs77a2h+vfuGXPZWHFH4+YOGbf15Wab3dR1tu4jbr4Xbr6Psaan8hiUDMk4I9gCDBWRIfgO4nOA/x0yzxrgbuB14IvARlVVEVkD/JeI/AjfxeKhwJtJiClI+wUV/90AIXcKBE6PdKdBpGnx9CcaadlNe46Fjbe9PNq627clVPt/H7G2Kdr09mnR7ghJxbLxfhaR6nrxzaPC3umy+OZRMT+n9vii3b0TSaz9K5pY7xsrbjfb7Ga/jzfudO8j8X4vIn2nwgn9PqbiWBDvZ5Ho/hUv1xeLAUTkJuCf8N0+ukxVHxORR4GtqrpGRHoBvwKKgJPAHFXd7yz7EPA14AzwLVV9Odb7+ZqGOnepoKKigqlTp3Zq2Wju+KXv+nYiv4IcsmAt4WpfgANL/hfguxMh3I4VeocC+P5LSPb96amqLzfu+OXrTt+tN3aYFuvWwc58TkHCNA0FysT+Fc82R6qvTEt2fcXzvYj0nYrn+xgP1/tYFG7rSyT8xeKkXCNQ1XXAupCyRwKGW4DbIyz7GPBYMuLoamI9eiCe2xKz9SfrmRL4H2eoVN+LnSrx/EahK2xHOsT6XkT7TmXro0DSoctcLO6O5s8YHva/l/a2v1jPe7EDQPzc3OufSV017kyK9r2I9p2K9X3szuwRExnU/uiB9vbL0EcPZOsD2Lqirvq4j64ad7aK9p2K9X3szuyMIMPa/3sJ1/bn5VPVZOuqSbWrxp2tYn2non0fuzM7I0iC9jbcvx44mdSHRXXHB3dlSld9GmZXjTtbpfo7lapjQapZInApUhtuMnaAwFNVwVunqsnWVZNqV407W6XyO5XKY0GqWdOQS2478IjFLggH6+ydP131Liu3cXfVO6VSKVXfqVQfC1LJEoFL1oabPm7voOmqSbWzcdsdR+nVlY8F1jTkkrXhpo/dQZMYq6/06srHAksELlkbbvp05f+4MsHqK7268rHAEoFLdkE3fbryf1yZYPWVXl35WGDXCJKgq7Y9dzVe/uVnZ1h9pV9XPRZYIjBdhpunfHqR1ZeJlyUC06V49ZefnWX1ZeJh1wiMMcbjLBEYY4zHWSIwxhiPc5UIROQCEXlFRPY5fweEmadQRF4XkV0islNE7giY9oyIHBCR7c6r0E08xhhjEuf2jGABsEFVhwIbnPFQnwBfUdVRwEzgn0Skf8D0+apa6Ly2u4zHGGNMgtwmglnAs87ws8Ds0BlU9R1V3ecMHwKOAhe7fF9jjDFJ4qrzehGpV9X+zrAAde3jEeYfjy9hjFLVsyLyDHAtcArnjEJVT0VYtgwoA8jLyytesWJFp2JuamqiT58+nVo2lSyuxGQqrqmlpQBUbNoUdrrVV2IsrsS4jau0tDRs5/WoatQX8CpQFeY1C6gPmbcuynoGAnuBa0LKBDgXX4J4JFY8qkpxcbF21qZNmzq9bCpZXInJWFzge0Vg9ZUYiysxbuMCtmqYY2rMH5Sp6vRI00TkiIgMVNXDIjIQX7NPuPk+BawFHlLVNwLWfdgZPCUi/wH8fax4jDHGJJfbawRrgLud4buBF0NnEJFc4AXgOVVdGTJtoPNX8F1fqHIZjzHGmAS5TQRLgOtFZB8w3RlHREpE5Glnnr8FpgBzw9wm+msRqQQqgYuA77uMxxhjTIJcPWtIVU8Anw9TvhW4xxn+T+A/Iyw/zc37G2OMcc9+WWyMMR5nicAYYzzOEoExxnicJQJjjPE4SwTGGONxlgiMMcbjLBEYY4zHWSIwxhiPs0RgjDEeZ4nAGGM8zhKBMcZ4nCUCY4zxOEsExhjjcZYIjDHG4ywRGGOMx7lKBCJygYi8IiL7nL8DIszXFtApzZqA8iEi8lcRqRGR3zi9mRljjEkjt2cEC4ANqjoU2OCMh9OsqoXO6+aA8ieAH6vqVUAdMM9lPMYYYxLkNhHMAp51hp/F1+9wXJx+iqcB7f0YJ7S8McaY5BBV7fzCIvWq2t8ZFqCufTxkvjPAduAMsERVy0XkIuAN52wAERkMvKyqoyO8VxlQBpCXl1e8YsWKTsXc1NREnz59OrVsKllciclUXFNLSwGo2LQp7HSrr8RYXIlxG1dpaelbqlrSYYKqRn0BrwJVYV6zgPqQeesirCPf+Xsl8B7wGXyd1dcEzDMYqIoVj6pSXFysnbVp06ZOL5tKFldiMhYX+F4RWH0lxuJKjNu4gK0a5pgas/N6VZ0eaZqIHBGRgap6WEQGAkcjrKPW+btfRCqAImAV0F9EzlHVM8AgoDZWPMYYY5LL7TWCNcDdzvDdwIuhM4jIABE51xm+CJgE7Hay0ybgi9GWN8YYk1puE8ES4HoR2QdMd8YRkRIRedqZZySwVUR24DvwL1HV3c60B4Fvi0gNcCHw7y7jMcYYk6CYTUPRqOoJ4PNhyrcC9zjDrwFjIiy/HxjvJgZjjDHu2C+LjTHG4ywRGGOMx1kiMMYYj7NEYIwxHmeJwBhjPM4SgTHGeJwlAmOM8ThLBMYY43GWCIwxxuMsERhjjMdZIjDGGI+zRGCMMR5nicAYYzzOEoExxnicJQJjjPE4SwTGGONxrhKBiFwgIq+IyD7n74Aw85SKyPaAV4uIzHamPSMiBwKmFbqJxxhjTOLcnhEsADao6lBggzMeRFU3qWqhqhYC04BPgD8EzDK/fbqqbncZjzHGmAS5TQSzgGed4WeB2THm/yLwsqp+4vJ9jTHGJImoaucXFqlX1f7OsAB17eMR5t8I/EhVf+eMPwNcC5zCOaNQ1VMRli0DygDy8vKKV6xY0amYm5qa6NOnT6eWTSWLKzGZimtqaSkAFZs2hZ1u9ZUYiysxbuMqLS19S1VLOkxQ1agv4FWgKsxrFlAfMm9dlPUMBI4BPUPKBDgX3xnFI7HiUVWKi4u1szZt2tTpZVPJ4kpMxuIC3ysCq6/EWFyJcRsXsFXDHFPPiZVBVHV6pGkickREBqrqYREZCByNsqq/BV5Q1daAdR92Bk+JyH8Afx8rHmOMMcnl9hrBGuBuZ/hu4MUo894JLA8scJJHe7PSbHxnGsYYY9LIbSJYAlwvIvuA6c44IlIiIk+3zyQiVwCDgT+GLP9rEakEKoGLgO+7jMcYY0yCYjYNRaOqJ4DPhynfCtwTMP4ekB9mvmlu3t8YY4x79stiY4zxOEsExhjjcZYIjDHG4ywRGGOMx1kiMMYYj7NEYIwxHmeJwBhjPM4SgTHGeJwlAmOM8ThLBMYY43GWCIwxxuMsERhjjMdZIjDGGI+zRGCMMR5nicAYYzzOEoExxnicq0QgIreLyC4ROSsiJVHmmykie0WkRkQWBJQPEZG/OuW/EZFcN/EYk0rl22qZ9HfLGPLAGiYt2Uj5ttpMh2RMUrg9I6gCbgX+FGkGEckBfgbcCBQAd4pIgTP5CeDHqnoVUAfMcxmPMSlRvq2Whasrqe33aVR6UFvfzMLVlZYMTLfgKhGoarWq7o0x23igRlX3q+ppYAUwy+mwfhqw0pnvWXwd2BuTdZau30tza1tQWXNrG0vXx9r9jcl+oqruVyJSAfy901dx6LQvAjNV9R5n/C5gArAYeMM5G0BEBgMvq+roCO9RBpQB5OXlFa9YsaJTsTY1NdGnT59OLZtKFldi0h3X3N9/HHHaMzPP9w9bfSXG4kqM27hKS0vfUtUOzfgxO68XkVeBS8JMekhVX+x0RAlS1aeApwBKSkp06tSpnVpPRUUFnV02lSyuxKQ7rvw3NlJb39yxvH/voDisvhJjcSUmVXHFTASqOt3le9QCgwPGBzllJ4D+InKOqp4JKDcm68yfMZyFqyuDmod698xh/ozhGYzKmORIx+2jW4Chzh1CucAcYI362qQ2AV905rsbSNsZhjGJmF2Uz+O3jiG/f28E35nA47eOYXZRfqZDM8a1mGcE0YjILcC/ABcDa0Vku6rOEJFLgadV9SZVPSMi9wHrgRxgmaruclbxILBCRL4PbAP+3U08xqTS7KJ8O/CbbslVIlDVF4AXwpQfAm4KGF8HrAsz3358dxUZY4zJEPtlsTHGeJwlAmOM8ThLBMYY43GWCIwxxuOS8svidBORY8B/d3Lxi4DjSQwnWSyuxFhcibG4EtNd47pcVS8OLeySicANEdka7ifWmWZxJcbiSozFlRivxWVNQ8YY43GWCIwxxuO8mAieynQAEVhcibG4EmNxJcZTcXnuGoExxphgXjwjMMYYE8ASgTHGeJxnEoGILBWRPSKyU0ReEJH+AdMWikiNiOwVkRlpjut2EdklImdFpCSg/AoRaRaR7c7rX7MhLmdaxuorJI7FIlIbUEc3xV4qpfHMdOqkRkQWZDKWQCLynohUOnXUoRfBNMaxTESOikhVQNkFIvKKiOxz/g7Ikrgyvm+JyGAR2SQiu53v4v9zypNfZ6rqiRdwA3COM/wE8IQzXADsAM4FhgDvAjlpjGskMByoAEoCyq8AqjJYX5Hiymh9hcS4GF8Xqdmwf+U4dXElkOvUUUGm43Jiew+4KAvimAKMC9yvgSeBBc7wgvbvZRbElfF9CxgIjHOG+wLvON+/pNeZZ84IVPUP6usJDeANfD2iAcwCVqjqKVU9ANSQxkdjq2q1qmZdD+hR4spofWWx8UCNqu5X1dPACnx1ZRyq+ifgZEjxLOBZZ/hZYHY6Y4KIcWWcqh5W1bed4UagGsgnBXXmmUQQ4mvAy85wPvBBwLSDTlk2GCIi20TkjyLyuUwH48i2+rrPae5blolmhQDZVi+BFPiDiLwlImWZDiZEnqoedoY/BPIyGUyIbNm3EJErgCLgr6Sgzlx1TJNtRORV4JIwkx5S1RedeR4CzgC/zqa4wjgMXKaqJ0SkGCgXkVGq2pDhuNIqWozAL4Dv4TvQfQ/4Ib4kb4JNVtVaEfk08IqI7HH+C84qqqoiki33s2fNviUifYBVwLdUtUFE/NOSVWfdKhGo6vRo00VkLvAF4PPqNLABtcDggNkGOWVpiyvCMqeAU87wWyLyLjAMSNrFvs7ERRrqK1C8MYrIvwG/S1UccUhrvSRCVWudv0dF5AV8zVjZkgiOiMhAVT0sIgOBo5kOCEBVj7QPZ3LfEpGe+JLAr1V1tVOc9DrzTNOQiMwEHgBuVtVPAiatAeaIyLkiMgQYCryZiRgDicjFIpLjDF+JL679mY0KyKL6cr4E7W4BqiLNmwZbgKEiMkREcoE5+Ooqo0TkfBHp2z6M76aJTNZTqDXA3c7w3UC2nIlmfN8S37/+/w5Uq+qPAiYlv84yeVU8zVfga/C14W53Xv8aMO0hfHd87AVuTHNct+BrTz4FHAHWO+W3AbucWN8G/iYb4sp0fYXE+CugEtjpfDkGZngfuwnfnR3v4mtey1gsATFdie8Oph3O/pSxuIDl+Jo8W519ax5wIbAB2Ae8ClyQJXFlfN8CJuNrmtoZcNy6KRV1Zo+YMMYYj/NM05AxxpjwLBEYY4zHWSIwxhiPs0RgjDEeZ4nAGGM8zhKBMcZ4nCUCY4zxuP8PfZxysnIf7Z8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "[xc,lags] = M_xcorr(x, y, maxlags=20)\n",
    "\n",
    "Xc = np.zeros(np.size(xc))\n",
    "Xc[19] = -1\n",
    "Xc[21] = 1\n",
    "\n",
    "plt.stem(lags,xc, label = '$Sample cross-correlation$')\n",
    "markerline, stemlines, baseline = plt.stem(lags,Xc, linefmt = 'r-', label = '$Theoretical autocorrelation$')\n",
    "plt.setp(stemlines, 'linewidth', 2)\n",
    "plt.grid(True)\n",
    "plt.legend(loc = \"upper left\")\n",
    "plt.axhline(0, color='blue', lw=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d04a4dc3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
